##test rbf solution

import RadialBasisFunctions as rbf
import numpy as np
import matplotlib.pyplot as plt
from math import log
from math import sqrt
from math import exp

s = 100
vol = 0.35
T = 3
r = 0.05
strike = 100
n = 51
p = 51


dt = T/float(p - 1)
xmax = log(s) + vol * sqrt(T) * 3.5
xmin = log(s) - vol * sqrt(T) * 3.5

x =  [ i * (xmax - xmin)/float(n-1) + xmin for i in range(n)]
dx = [x[i] - x[i-1] for i in range(1,len(x))]
c = 4* min(dx)

v = np.array([max(exp(x[i]) -strike,0.0) for i in range(len(x)) ])

f = lambda x,xc,m : rbf.MQ([x],[xc],c,m)
##initialise coeffs
A = np.array([[f(x[i],x[j],0) for  j in range(len(x))] for i in range(len(x))])

B = np.array([[f(x[i],x[j],1) for  j in range(len(x))] for i in range(len(x))])

C = np.array([[f(x[i],x[j],2) for  j in range(len(x))] for i in range(len(x))])

##matrix to solve
D  = 0.5* vol* vol* C  + ( r - 0.5*vol* vol)* B - (r + 1.0/dt) * A

for k in range(1,p-1):
    coefs = np.linalg.solve(D,-v/dt)
    ##update v
    v = np.array( [ np.sum(coefs[j] * f(x[i],x[j],0) for j in range(len(x))) for i in range(len(x)) ])
    v[0] = 0
    v[-1] = exp(x[-1]) -strike
    coefs = np.linalg.solve(A,v)


print "Price at spot "  + str( np.sum([f(log(s),x[i],0)*coefs[i] for i  in range(len(x))]) )
delta  = np.array( [ np.sum(coefs[j] * f(x[i],x[j],1)/exp(x[i]) for j in range(len(x))) for i in range(len(x)) ])
gamma  = np.array( [ np.sum(coefs[j] * f(x[i],x[j],2)/exp(2*x[i]) for j in range(len(x))) for i in range(len(x)) ])

fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_ylabel('price')
ax.set_xlabel('log S')
ax.plot(x,v)

ax2 = ax.twinx()
ax2.plot(x,delta)
ax2.set_ylabel('delta')
plt.show()





